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1.  Summary 


Multidimensional  numerical  simulation  (three-dimensional-computational  fluid  dynamics  [3D- 
CFD])  of  the  spray  behavior  is  an  integral  part  of  engine  research  critical  for  combustion 
optimization  and  to  enhance  the  understanding  of  complex  in-cylinder  combustion  phenomena 
occurring  in  engines.  Successful  implementation  of  advanced  modeling  tools  is  strongly  related 
to  our  current  understanding  of  the  physical  processes  involved.  The  importance  of  numerical 
models  and  simulation  tools  is  increasing  dramatically  (7).  Recently,  it  has  been  essential  in 
aiding  the  development  process  of  new  combustion  concepts  providing  unprecedented 
efficiencies  (2).  Such  developments  are  critical  to  engine  technology  advancements  of  interest  to 
the  U.S.  Anny.  The  objective  of  this  report  is  to  review  the  leading  3D-CFD  atomization  models 
used  in  engine  research  to  illustrate  our  current  understanding,  and  to  identify  grey  areas  of 
further  research. 

Spray  behavior  can  be  described  as  a  multiscale,  turbulent  physical  process  presenting  several 
technical  challenges.  There  is  a  liquid  core  region  that  is  affected  by  the  aerodynamic  interaction. 
Once  the  liquid  surface  becomes  unstable  it  will  favor  the  creation  of  liquid  ligaments  that  in  turn 
will  create  parent  “primary”  droplets  at  first,  followed  by  secondary  “child”  droplets.  Droplets 
are  reduced  in  size  due  to  evaporation  and  combustion  occurs  while  reduced  droplets  traveling 
downstream  of  the  injector  nozzle.  Models  capturing  this  process  have  had  success  in  engineering 
applications,  and  have  typically  relied  on  empirical  correlations  or  simplified  physics.  However, 
deeper  knowledge  about  turbulent  sprays  is  still  needed  in  the  community.  Several  reports  have 
noted  an  insufficient  understanding  of  the  development  of  surface  instabilities,  ligament  fonnation 
mechanisms,  and  droplet  formation  mechanism  from  a  ligament.  Nozzle  effects  and  needle 
wobble  conditions  are  also  important  characteristics  that  have  not  been  fully  resolved  and 
strongly  affect  spray  breakup.  Improved  knowledge  of  primary  breakup  will  lead  to  better 
predictions  of  spray  characteristics  such  as  initial  droplet  size  distribution,  spray  angle,  and  jet 
structure,  which  result  in  significant  impacts  on  improving  the  engine  combustion  process. 

Spray  breakup  models  may  be  improved  through  the  use  of  high-resolution  numerical 
simulations.  Large  Eddy  Simulation  (LES)  can  provide  a  tradeoff  between  accuracy  and 
computational  cost.  LES  is  dependent  on  the  subgrid  models  to  treat  unresolved  flow-field. 

Direct  Numerical  Simulation  (DNS)  is  aimed  at  removing  turbulence  uncertainty  and  resolving 
all  of  the  flow  length-scales.  Of  key  importance  is  resolving  the  droplet  length  scale,  in  addition  to 
the  flow  length  scale,  and  studying  the  droplet  interactions  with  its  surroundings.  Computationally, 
DNS  is  still  a  daunting  task  due  to  its  resolution  requirements  imposed  by  the  droplet  scale. 
Recently,  several  studies  using  DNS  tools  have  made  significant  strides  in  primary  atomization 
models  using  the  Eularian  Lagrangian  Spray  Atomization  (ELSA)  methodology  and  driving  the 
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development  of  advanced  models.  Typically,  these  studies  are  conducted  in  small  configurations, 
and  are  representative  of  efforts  towards  next  generation  atomization  models  for  diesel  sprays. 


2.  Introduction 


The  objective  of  this  report  is  to  review  contemporary  liquid  spray  physical  models  used  in  the 
computational  fluid  dynamics  (CFD)  community  for  engine  applications.  It  provides  a  technical 
introduction,  a  presentation  of  the  physical  mechanisms,  and  a  review/discussion  of  the  leading 
models  and  applications. 

The  aims  of  the  report  are  as  follows: 

•  To  provide  a  review  of  current  spray  modeling  efforts  used  in  the  state-of-the-art  engine 
research. 

•  To  identify  critical  areas  of  future  modeling  interest  that  can  lead  to  technological 
advancements. 

•  To  promote  engine  spray  research  by  dissemination  of  these  findings  to  the  scientific 
community. 

Fundamental  physical  understanding  of  the  spray  atomization  process  is  (1)  critical  in  addressing 
the  Army’s  need  to  enable  the  advancement  of  efficient  air  and  ground  combat  vehicles  and 
compact  fuel-to-electric  power  systems,  and  (2)  essential  to  continue  leading  strategic 
development  for  portable  and  heavy-duty  engine  applications.  Spray  technology  has  a  wide 
range  of  applications  in  engineering — such  as  combustion  devices  found  in  unmanned  aerial 
systems  (rotary  and  piston  engines),  aerospace  propulsion  systems  (gas  turbine  engines),  ground 
vehicle  engines  (piston  engines),  metallurgical  processes  (spray  generated  coatings)  for  thermal 
barrier,  environmental  barrier,  substrate  corrosion  protection,  and  material  strengthening  (3-5). 
Within  the  context  of  unmanned  aerial  and  ground  vehicle  engine  applications,  an  understanding 
of  liquid  spray  characteristics  is  critical  in  improving  engine  perfonnance  and  efficiency  as  well 
as  complying  with  strict  emissions  standards  that  require  precise  fuel  delivery  systems  for 
optimum  power  and  reduction  of  exhaust  gases  and  particulates.  Hence,  detailed  characterization 
of  the  spray  and  combustion  processes  is  important  in  developing  high-impact  new  technologies 
such  as  single-fuel  optimized  engines  for  Army  engine  applications. 

Spray  breakup  is  a  complex  phenomenon  and  dependent  on  many  factors  such  as  velocity  of 
spray,  injection  pressure,  geometry  of  nozzle,  temperature  of  both  fluids,  and  properties  of  fluid 
(i.e.,  their  density,  viscosity,  surface  tension,  and  vapor  pressure)  (6).  Several  detailed  studies  of 
the  spray  process  exist  in  the  combustion  literature  focusing  on  the  atomization  and  breakup 
characteristics  of  fuel  sprays  through  the  use  of  sophisticated  combustion  vessels  (7-10).  The 
Engine  Combustion  Network  (ECN)  established  through  the  efforts  at  Sandia  National 


2 


Laboratories  (SNL)  has  historically  driven  the  experimental  efforts  in  quantifying  spray 
parameters  with  the  aim  of  providing  a  high-fidelity  consolidated  database  at  diesel  and  gasoline 
engine  relevant  conditions.  Some  of  the  ECN  affiliates  that  have  conducted  spray 
characterization  in  vessels  include:  SNL,  French  Petroleum  Institute  (IFPEN),  Universidad 
Politecnica  de  Valencia  (CMT-Motores  Termicos),  Caterpillar,  Eindhoven  University  of 
Technology  (TU/e),  Michigan  Technological  University  (MTU),  and  General  Motors.  Most 
recently,  the  U.S.  Army  Research  Laboratory  (ARL)  is  expanding  its  involvement  in  the  ECN 
network.  Measurements  are  typically  conducted  at  well  controlled  high  temperature  and  pressure 
conditions  through  the  use  of  two  distinct  types  of  combustion  vessels:  Constant-Pressure  Flow 
(CPF)  and  Constant- Volume  Prebum  (CVP)  chambers.  An  excellent  review  has  been  published 
focusing  on  the  comparison  between  the  operations  and  control  of  boundary  conditions  for  both 
types  of  combustion  vessels  (CPF  and  CVP)  while  studying  a  well  reported  case  denoted  as 
Spray  A  (11).  Spray  A  uses  a  single  component  diesel  surrogate  fuel  (n-dodecane),  a  single  hole 
injector  (common  rail,  1500 -bar  fuel  pressure ,  363 -K fuel  temperature),  representing  a  diesel 
engine  combustion  condition  (900  K,  60  bar)  that  use  a  moderate  rate  of  exhaust-gas 
recirculation  (EGR).  Using  a  CVP  chamber,  Siebers  et  al.  (12)  reported  the  spray  liquid 
penetration  length  at  near  Spray  A  conditions  and  studied  the  spray  behavior  when  subject  to 
changes  in  operating  conditions  such  as  decreasing  injector  orifice  diameter,  injection  pressure, 
ambient  gas  density  or  temperature,  and  fuel  volatility.  It  was  established  that  liquid  length 
decreases  linearly  with  injector  diameter,  temperature  or  density,  increases  with  fuel  volatility  or 
temperature,  and  that  it  is  independent  of  injection  pressure.  Extending  this  work,  Weber  et  al. 
(13)  has  used  a  CPF  type  vessel  at  diesel  conditions,  together  with  computational  methods  (CFD 
and  microgenetic  algorithms),  to  provide  optimization  strategies  for  spray  penetration  and 
mixture  formation  at  diesel  condition  (50  bar  and  800  K).  Similar  research  efforts  led  by  Kweon 
(3)  have  been  performed  at  the  ARL  for  diesel  and  JP-8  fuels  to  characterize  the  effects  of 
injector  configuration  and  fuel  composition  on  the  spray  parameters  using  a  CPF  vessel  and 
optical  diagnostics  such  as  Mie  scattering,  Shadowgraph,  and  Schlieren  images. 

Access  to  spray  measurements  plays  a  key  role  in  the  development  of  physics-based  models  to 
be  integrated  into  robust  3D-CFD  solvers.  Multidimensional  CFD  solvers  will  utilize 
experimental  data  by  carrying  out  model-validation  studies  and  quantifying  the  simulation  error; 
hence  the  solver  strictly  relies  on  the  fidelity  of  the  measurement  in  addition  to  its  own  numerical 
errors.  A  validation  study  will  reveal  the  suitability  of  modeling  assumptions  (physical  models), 
stability  of  the  spatio-temporal  numerical  technique  (numerical  methods),  and  calibration  of 
model  parameters  (turbulence,  breakup,  combustion  constants)  that  is  required  to  carefully 
benchmark  a  simulation.  Hence,  validated  modeling  of  spray  and  combustion  processes  enables 
several  attributes  that  make  its  use  in  engine  development  of  paramount  importance;  for 
example,  its  suitability  in  conducting  extensive  cost-effective  parametric  studies,  the  ability  to 
provide  a  complete  history  of  multidimensional  resolved  infonnation  for  every  variable,  and  its 
ability  to  artificially  separate  specific  subprocesses  that  would  inherently  interact  in  an 
experiment  to  provide  new  insights.  Over  the  years,  these  attributes  have  contributed  to  the  use 
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of  CFD  solvers  in  spray  processes  for  both  research  and  design  purposes.  Some  of  the  CFD 
solvers  representing  the  current  state  of  the  art  of  spray  modeling  are:  CONVERGE  developed 
by  Convergent  Science,  USA  (14);  KIVA,  developed  by  Los  Alamos  National  Laboratories, 

USA  (15);  OpenFoam  developed  by  OpenCFD,  U.K.;  and  AVBP  developed  by  Centre  Europeen 
de  Recherche  et  de  Formation  Avancee  en  Calcul  Scientifique  (CERFACS),  France.  The  solvers 
are  available  for  the  simulation  of  complex  combustion  phenomena  and  are  representative  of 
ongoing  efforts  within  the  combustion  community  to  advance  the  general  area  of  CFD-based 
engine  research. 

In  recent  years,  the  engine  modeling  community  has  utilized  SNL  ECN  network  spray  data 
mainly  to  validate  and  demonstrate  several  numerical  features  of  CFD  solvers  (16-19).  These 
features  include  the  solvers  ability  to  model  turbulence  in  the  flow  field,  obtain  grid-convergent 
simulations,  and  study  the  effect  of  ensemble  averaging  on  spray  kinematics,  penetration  length, 
and  other  important  parameters  (14, 17).  Modeling  the  transients  in  the  flow  field  is  an  important 
feature  due  to  the  spray  local  inhomegeneities  that  may  arise  during  mixing.  Traditionally,  coarse 
turbulence  models  have  been  used  in  engine  research  due  to  the  burden  of  the  computational  cost 
associated  with  the  grid-resolution.  Turbulence  modeling  is  classified  on  its  level  of  flow/grid 
resolution  and  its  cost,  with  DNS  resolving  all  the  length  scales,  LES  resolving  the  anisotropic 
lengthscales  while  modeling  the  isotropic/dissipation  scales,  and  Reynolds  Average  Navier 
Stokes  (RANS)  that  is  based  on  ensemble  averaging  of  the  governing  equations.  Using  a  RANS 
technique,  Senecal  et  al.  (1 7)  has  reported  in  detail  the  modeling  guidelines  to  simulate 
nonevaporating,  evaporating,  and  reacting  sprays  utilizing  the  CONVERGE  solver.  In  this  study, 
ECN  spray-A  was  used  to  demonstrate  key  modeling  features  such  as  the  importance  of 
adaptive-mesh-refinement  (AMR),  grid  embedding,  grid  convergence  properties,  and  critical 
lagrangian  parcel  distribution  for  accurate  spray  modeling  (17).  Further  studies  demonstrated  the 
effect  of  several  high-resolution  turbulence  methods  (LES  techniques)  on  resolving  the  flow 
field,  cycle-to-cycle  variations  in  the  context  of  grid  convergence,  and  its  effect  on  spray 
modeling  constants  (18,  19).  Good  agreement  was  found  for  global  quantities  such  as  liquid  and 
vapor  penetration  lengths  with  experimental  data  for  single-shot  LES  realization.  However, 
deeper  comparison  to  local  quantities  such  as  velocity  and  mixture  fraction  was  shown  to  require 
the  average  of  many  LES  realizations.  Additionally,  the  spray  breakup  constants  are 
demonstrated  to  depend  on  the  amount  of  resolution  utilized  in  resolving  the  flow  fields,  and  the 
transition  from  RANS  to  LES  turbulence  methods  does  require  consideration  to  recalibrate. 
Similarly,  Habchi  and  Bruneaux  (20)  studied  the  effects  of  LES  ensemble  averaging  for  a  single 
hole  injector  in  a  high-pressure,  high-temperature  vessel  and  found  good  agreements  with 
experimental  data  using  10-30  LES  realizations  using  the  AVBP  solver  (20).  Hence,  it  is 
generally  agreed  that  using  LES  is  similar  to  a  single-shot  experimental  injection  requiring 
ensemble  averaging,  while  with  RANS  a  single  shot  realization  will  be  enough  to  capture 
important  spray  quantities. 
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In  the  preceding  text,  several  state-of-the-art  methodologies  have  been  presented  highlighting 
recent  progress  in  both  experimental  and  computational  engine  research.  In  the  context  of  spray 
modeling,  the  predictive  quality  of  the  phenomelogical  models  in  CFD  solvers  have  reached  a 
mature  level  making  its  use  a  necessity  in  engine  research  and  development  (7,  2,  21).  In  spite  of 
such  progress,  there  are  several  areas  that  have  not  been  fully  addressed  by  the  modeling 
community  and  where  progress  is  clearly  needed.  Phenomelogical  breakup  models  that 
accurately  account  for  primary  breakup  effects  (near  the  nozzle  orifice)  have  not  been  fully  able 
to  simulate  the  process.  It  is  also  known  that  spray  is  optically  very  dense  in  this  region,  making 
it  difficult  to  measure.  As  a  result,  a  comprehensive  understanding  of  the  contributions  of 
cavitation,  multiphase  turbulence,  and  aerodynamics  effects  on  the  spray  process  has  not  been 
achieved  presently.  Additionally,  criteria  for  grid-converged  spray  calculations  are  still  being 
studied  and  no  globally  accepted  consensus  exists  today;  however,  much  progress  has  been 
reported  for  high-resolution  models. 

The  reports  cover  the  following  areas: 

•  An  introduction  to  existing  methodologies  (both  experimental  and  computational)  used  to 
carry  out  fundamental  spray  research  in  engine  applications. 

•  A  review  of  key  physical  mechanisms,  and  governing  parameters,  of  importance  in  spray 
model  development  and  analysis. 

•  A  review  of  contemporary  liquid  spray  models  used  in  CFD  modeling  community  and 
provide  model  application  guidelines. 

•  A  presentation  of  technical  challenges  and  unexplored  areas  in  spray  research  providing 
outlooks  and  recommendations  in  engine  spray  modeling. 


3.  Physical  Description 

Sprays  play  the  role  in  mixing  the  fuel  with  air,  and  to  increase  its  surface  area  for  rapid 
evaporation  and  combustion.  Its  process  significantly  affects  the  ignition  behavior,  heat  release 
and  pollutant  formation  rates,  fuel  consumption,  and  exhaust  emissions.  The  kinetic  energy  of 
the  spray  represents  the  main  source  for  turbulence  production  and  governs  the  microscale  air- 
fuel  mixing  by  turbulent  diffusion  and  the  flame  speed  of  the  premixed  flame  front  (21). 

The  spray  and  atomization  processes,  can  be  described  as  a  multiphase  flow  phenomena 
involving  liquid  phase  in  the  form  of  droplets  and  ligaments  and  the  gas  phase  represented  as  a 
continuum.  The  spray  process  is  typically  initiated  as  a  result  of  high-pressure  liquid  fuel 
discharged  from  an  injector  nozzle  (or  several  nozzles)  carrying  with  its  important  physical 
features  like  liquid-phase  turbulent  flows,  and  cavitation  effects  from  the  generation  of  gas-phase 
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bubbles  that  can  potentially  implode  as  they  travel  internally  downstream  of  the  flow  and  are 
ejected  into  the  chamber.  With  reference  to  the  above  discussion,  figure  1,  provides  an 
appropriate  description  of  the  diesel  spray  atomization  process. 


Figure  1.  Various  stages  of  high-pressure  diesel  spray  breakup. 

The  disintegration  or  breakup  occurs  when  the  disruptive  forces  exceed  surface  tension  forces. 
External  forces — such  as  aerodynamic  forces,  surface  shear  forces,  centrifugal  forces  and 
electrostatic  forces,  acting  on  the  liquid  surface  may  distort  the  bulk  liquid  and  promote  the 
disruption.  External  forces  may  lead  to  oscillations  and  perturbations  of  the  interfaces.  These 
oscillations  may  be  amplified  and  result  in  the  breakup  of  the  liquid  into  small  droplets.  This 
initial  breakup  process  is  often  referred  to  as  the  primary  breakup  or  the  primary  atomization.  A 
population  of  larger  droplets  produced  in  the  primary  atomization  may  be  unstable  if  they  exceed 
a  critical  droplet  size  and  thus  may  undergo  further  disruption  into  smaller  droplets.  This  process 
is  usually  termed  as  the  secondary  breakup  or  the  secondary  atomization.  Therefore,  the  final 
droplet  size  distribution  produced  in  an  atomization  process  is  detennined  by  the  flow 
characteristics  and  the  properties  of  the  fluids  in  both  the  primary  and  secondary  disintegration. 

If  the  surrounding  temperature  is  high  enough,  the  droplets  will  evaporate  producing  vapor 
which  mixes  with  the  oxidizer  forming  a  combustible  mixture  which  ignites  due  to  the  presence 
of  sparks,  or  due  to  increased  pressures  and/or  temperature  in  compression-ignition  engines. 

The  analysis  of  atomization  and  sprays  are  typically  carried  out  by  means  of  theoretical, 
numerical,  or  experimental  methodologies  (22).  As  in  traditional  fluid  mechanics,  the 
characterization  of  spray  behavior  is  also  most  conveniently  analyzed  with  several  non¬ 
dimensionless  parameters. 
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Several  of  the  important  parameters  are  listed  below. 

Liquid  Reynolds  number.  pc  _Piuidi 

1  Pi 

(1) 

Liquid  Weber  number. 

...  Piufdt 

Wel  = - 

a 

(2) 

Aerodynamic  (gas)  Weber 
number. 

=  pXeA 
y  o 

(3) 

Ohnesorge  number. 

ohc 

\wei  _  Pi 

(4) 

Rei  y/p^di 

Discharge  Coefficient : 

fa 

vl 

12  AP 
i  Pi 

(5) 

Cavitation  Parameter: 

K  = 

.2  (Pl-Pg) 

PlV 2 

(6) 

3.1  Primary  Atomization 

In  primary  breakup,  the  behavior  of  liquid  sheets  (or  jets)  can  be  classified  into  different 
atomization  modes  depending  on  boundary  conditions.  Figure  2  depicts  the  droplet  breakup 
behavior  with  increasing  Weber  number  conditions.  Earlier  investigation  of  single  fluid  pressure 
atomization,  divided  the  breakup  regimes  of  a  circular  liquid  jet  into  three  areas  depending  on  the 
liquid  Reynolds  number  and  the  Ohnesorge  number  (23).  The  regimes  are  further  described 
below,  and  are  also  illustrated  in  figure  3. 

1.  At  low-Reynolds  number,  the  jet  disintegrates  due  to  surface  tension  effects  into  fairly 
identical  droplet  sizes  (Rayleigh  regime,  symmetric,  or  varicose  instability). 

2.  At  intermediate-Reynolds  number,  drop  formation  is  influenced  by  aerodynamics  forces 
(non-axisymmetric  Rayleigh  breakup).  These  forces  cause  symmetric  (first  wind-induced 
mode)  and  asymmetric  (second  wind-induced  mode,  asymmetric  or  sinusoidal  instability) 
wave  growth  of  gas  liquid  interface  that  finally  leads  to  jet  disintegrations.  This  is  known 
as  the  aerodynamic  regime). 

3.  At  higher-Reynolds  number,  the  jet  disintegrates  almost  spontaneously  at  the  nozzle  exit. 
This  is  called  the  atomization  regime. 


The  transition  line  from  the  Rayleigh  regime  to  the  aerodynamic  mode  is  identified  by 


Wel  = 


1.74  x  104 

Re °-5 


(V) 
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and  the  transition  line  between  aerodynamic  breakup  and  the  atomization  regime  is  characterized 
as 


Wel 


9.4  x  10s 
Re0-5 


(8) 


Figure  2.  Disintegration  of  circular  water  jet  in  pressure  atomization  (d=1.5  mm)  (24). 


Figure  3.  Primary  fragmentation  modes  of  a  liquid  jet  in  pressurized  atomization  (23). 
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The  atomization  process  is  also  characterized  by  the  continuous  jet  breakup  length,  where  the 
fuel  remains  as  a  continuous  medium.  It  is  the  distance  from  the  nozzle  exit  to  the  breakup  point, 
the  general  behavior  of  the  breakup  length  and  its  dependence  on  the  jet  velocity  as  shown  in 
figure  4. 


Figure  4.  Jet-stability  curve  showing  the  breakup  length  modes  and  dependence  with  velocity. 

The  initial  part  of  the  curve  is  described  as  the  dripping  region  of  the  jet.  The  laminar  flow 
region  is  located  from  point  A  to  point  B,  where  symmetric  Rayleigh  instabilities  prevail.  The 
upper  point  B  indicates  the  transition  from  varicose  to  sinusoidal  breakup  mode.  The  breakup 
length  decreases  in  the  transition  region  B  to  C.  When  the  fluid  at  the  nozzle  exit  is  already  in  a 
turbulent  flow  stage,  and  aerodynamic  interaction  between  the  liquid  jet  and  the  gas  dominates 
the  breakup,  the  jet  breakup  length  increases  with  increasing  velocity  (from  point  C  to  D).  The 
behavior  of  the  liquid  jet  breakup  length  at  jet  velocities  beyond  point  D  is  not  uniquely  defined 
yet,  but,  in  general,  tends  to  decrease.  Several  breakup  length  equations  have  been  proposed 
empirically  and  are  provided  in  textbook  references  {24). 

3.2  Secondary  Atomization 

In  secondary  atomization,  the  behavior  is  defined  as  the  disintegration  of  larger  droplets  and 
ligaments  into  smaller  droplets.  The  breakup  in  a  single  droplet  can  be  caused  by  a  relative 
velocity,  turbulence,  heat  and/or  mass  transfer.  Moreover,  the  instabilities  that  occur  for  high 
relative  velocities  between  the  deformable  liquid  droplet  and  the  surrounding  fluid  mainly  cause 
secondary  fragmentation  (22).  For  very  low-relative  velocities  the  droplet  remains  stable.  If  the 
aerodynamic  forces  overcome  the  forces  due  to  surface  tension,  the  droplet  will  tend  to  deform. 

For  higher  relative  velocities  the  deformed  droplet  will  breakup.  The  gas  or  aerodynamic  Weber 
number  shown  in  equation  3  (W eg)  gives  the  main  criterion  for  breakup.  If  Weber  numbers 
exceed  a  critical  Weber  number  Wecri  the  droplet  breaks  up.  As  surface  tension  has  a  stabilizing 
effect,  an  increase  in  viscosity  damps  unstable  perturbations.  The  effect  of  viscosity  on  the 
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breakup  is  included  in  the  Ohnesorge  number  as  in  equation  4.  For  Ohnesorge  numbers 
Oh  <  0.1,  the  effect  of  viscosity  is  negligible.  According  to  previous  reviews  on  droplet,  the 
behavior  of  the  Weber  number  has  been  shown  to  follow, 

(constant  for  Oh  <  0.1  and 

^ 6  l  ~Oh 2  for  Oh  >  0.1  (9) 

The  critical  Weber  number  found  for  low-viscosity  liquids  ranges  from  6  to  13,  and  it  is  usually 
reported  as 


Wecr  =  12 


(10) 


Viscosity  effects  increase  for  Weber  numbers  above  Wecr  and  different  types  of  droplet  regimes 
are  obeyed.  The  classification  of  breakup  regimes  includes  a  description  of  breakup  mechanisms 
for  high  Weber  number  Wei  —  350. 

Two  types  of  breakup  mechanisms  in  sprays  have  been  commonly  reported  as:  the  bag  breakup 
for  low-Weber  numbers  and  the  shear  breakup  for  high  Weber  numbers  as  shown  in  figure  5. 
Bag  breakup  is  typically  associated  with  Kelvin  Helmholtz  (KH)  instability  (parallel  shear  flow), 
while  the  shear  breakup  at  high  Weber  numbers  is  associated  with  the  Rayleigh  Taylor  (RT) 
instability  (cross  flow).  The  dynamics  of  secondary  fragmentation  process  depends  on 
aerodynamic  viscous  forces  (25).  For  spray  analysis,  in  a  simplified  model  of  the  dynamic 
breakup  process,  the  droplet  for  impulsive  aerodynamic  loading  initially  deforms  to  a  disc  and 
then  breaks  up  into  new  droplets.  A  characteristic  time  scale  needed  for  breakup  is  given  as, 


d0  Pi 
urel  Pg 

N 


(11) 
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High  Weg :  Stripping  breakup 


Figure  5.  Breakup  modes  in  secondary  droplet  breakup. 

For  computations  of  the  droplet  breakup  process,  CFD  spray  models  will  include  the  following 
criteria  for  determining  droplet  behavior. 

•  Criteria  to  determine  if  the  droplet  breaks  up  or  not 

•  An  estimation  of  the  breakup  time  tbup 

•  The  resulting  droplet  size  distribution  d32 

•  The  spatial  distribution  of  the  resulting  droplets  and  their  velocities  after  breakup 


4.  Spray  Modeling 


In  spray  modeling,  two  methodologies  that  have  received  the  most  application  are  the  Eularian 
approach,  and  the  Lagrangian  approach  (26,  27).  In  the  Eularian  approach,  the  liquid  is  treated  as 
a  continuous  medium  or  as  the  second  phase  of  the  multiphase  flow.  In  Lagrangian 
implementation,  the  modeling  is  focused  on  tracing  the  trajectories  of  individual  droplets  or 
parcel  of  droplets  alongside  the  continuum.  One  of  the  salient  features  of  Lagrangian  spray 
models  is  their  focus  on  jet  and  droplet  breakup.  Disintegration  of  the  liquid  jet  occurs  near  the 
nozzle  exit  in  the  primary  breakup  region.  Far  downstream  of  the  nozzle,  secondary  breakup  of 
large  droplets  into  smaller  ones  takes  place.  As  was  discussed  in  earlier  sections,  in  these  models 
the  jet  is  approximated  by  a  chain  of  droplets,  with  initial  diameters  equal  to  the  diameter  of  the 
nozzle,  or  slightly  less  than  this  diameter  if  the  effects  of  cavitation  are  taken  into  account. 


11 


This  section  aims  at  providing  several  models  for  separate  treatment  of  these  two  regions  or 
models  for  combined  treatment. 

4.1  The  Spray  Equation:  A  Statistical  Description 

A  statistical  description  for  the  stochastic  behavior  of  spray  was  introduced  on  the  landmark 
paper  by  Williams  (28).  This  mathematical  model  includes  the  effects  of  droplet  growth,  the 
formation  of  new  droplets,  collisions,  and  aerodynamics  forces.  In  principle,  for  sprays 
composed  of  N  particles,  a  statistical  description  will  involve  the  specification  of  the  Nth  order 
distribution  function  in  the  phase  space  consisting  of  3  A  spatial  coordinates  and  3N  velocity 
coordinates.  In  the  probability  density  approach,  the  spray  droplets  are  represented  by  a 
PDF,/(t,  x,  y,  z ),  representing  the  probable  number  of  droplets  per  unit  volume  in  time  and 
space.  The  state  of  the  droplet  is  described  by  its  parameters  that  are  the  coordinates  in  the 
particle  state  space,  the  parameters  typically  are:  location  x,  velocity  v,  radius  r  ,  the 
temperature  Td  ,  the  defonnation  parameter  and  rate  of  defonnation,  y,  y,  respectively. 

The  spray  PDF  is  the  solution  of  a  spray  transport  equation  given  in  component  form  as, 
df  d  d  d  d 

—  +  divx(fv)  +  clivy(fv)  +  —  (ff)  +  —  (ftd)  +  —  (fy)  +  —  (fy)  =  fcM  +  fbu  (12) 


where  the  source  terms  are  due  to  droplet  collision  and  breakup  and  are  modeled  via  individual 
spray  sub  processes.  Additional  models  also  need  to  be  available  for  each  subprocess  on  the  left 
hand  side  of  equation  10  for  a  complete  description  of  the  spray  physics. 

In  this  next  section,  models  are  presented  to  account  for  droplet  acceleration,  droplet  radius, 
droplet  temperature,  droplet  defonnation,  droplet  collision,  and  droplet  breakup  (atomization)  for 
use  in  contemporary  spray  solvers  (29). 


4.1.1  Droplet  Acceleration 


The  droplet  acceleration  equation  of  motion  is  given  by, 


— v  =  -C, 
dt 


Pg  llVrll 
8  r 


vr  +  g 


(13) 


where  Cd  is  the  drag  coefficient,  pg,  pL  the  gas  and  liquid  densities,  vr  is  the  relative  drop 
velocity  defined  as  vr  =  u  +  u'  —  v  accounting  for  turbulence  fluctuations,  and  g  is  the 
gravitational  constant.  The  drag  coefficient  is  given  by  following  formulas, 

c*4i^(1  +  Re‘2/3/6)  lf  Re“ £  1000 

(0.424  if  Red  >  1000 

where  Red  is  the  droplet  Reynolds  number  defined  by,  Red  —  2 rpgvr/pg(T),  and  the  weighted 
gas  temperature  is  given  by  the  two  third  law,  T  =  (T  +  2Td)/3. 
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4.1.2  Drop  Radius 


Considering  mainly  for  convective  heat  and  mass  transfer,  the  Frossling  correlation  is  used  to 
detennine  the  droplet  radius.  It  is  detennined  from  the  mass  rate  of  change  due  to  evaporation  or 
condensation,  its  expression  is  written  as, 

"wr2  =  —  DvBdShd  (14) 

at  pd 

where  Dv  is  the  vapor  diffusivity  in  the  gas  and  it  is  detennined  from  the  empirical  Frossling 
conelation,  pvDv  —  D^f0'1  having  D1  and  Dz  as  constants  and  T.  The  Spalding  mass  transfer 
number  is  used  to  define, 

Bd  =  (Y;  -  Yv)/(l  -  Yp)  (  . 


and  Yv  —  pv/pg  is  the  vapor  mass  fraction,  and  Fv*  is  the  vapor  mass  fraction  on  the  drop  surface 
calculated  assuming  equilibrium  conditions  and  invoking  the  Clayperon  thennodynamic 
equation, 


YvO'd) 


1  + 


MW,/  Pg 
MWV  \pv(Td ) 


(16) 


The  molecular  weights  are  denoted  as  MWS,  for  the  surrounding  gas,  and  MWV  for  the  vapor.  The 
equilibrium  vapor  pressure  is  denoted  as,  pv(Td)  and  pg is  the  gas  pressure. 


The  Sherwood  number  is  denoted  as 

Shd  =  (2.0  +  0.6 Re]/2 Sc]/3  )ln1  +  Bd 

Bd 


(17) 


where  the  droplet  Schmidt  number  is  defined  as,  Scd  —  p(f  )/ pgDg(f  ). 

4.1.3  Drop  Temperature 

An  energy  balance  determines  the  rate  of  change  of  the  droplet  temperature,  where  the  energy 
supplied  to  the  drop  raises  the  temperature  or  supplies  heat  for  evaporation.  The  mathematical 
expression  is  written  as, 

dT d  cLmd  . , 

Cdmd~dt  =qhSd+L(-Td)~dt  (18) 

where  Cd  is  the  droplet  specific  heat,  md  the  mass,  Sd  the  surface,  qh  is  the  convective  heat 
transferred  to  the  drop,  and  L(Td)  is  the  latent  heat  of  evaporation.  The  change  in  droplet 
temperature  heat  conduction  rate,  qh  ,  is  given  by  the  Ranz-Marshall  correlation  as  follows, 
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qh  =  Kg(T)Nud(T-Tdy2r 


(19) 


The  conductivity  of  the  gas  is  given  by  the  empirical  relationship  Kg(T )  =  ^ K \T  2J/(f  +  K2 ), 

where  Kx  and  Kz  are  constants.  The  convective  heat  transfer  is  expressed  in  terms  of  Nusselt 
number  as, 

Nud  =  (2.0  +  0.6Red/2Pr1d/3  )ln1  +  Bd  (20) 

&  d 

where  Red is  the  droplets  Reynolds  number  and  the  droplet  Prandlt  number  is  defined  by 
Prd  =  fi  g  (T )  Cp  (T)  / Kg  (f ) .  Finally  the  latent  heat  of  vaporization  is  defined  as  the  difference 
between  the  vapor  enthalpy  hv  and  the  liquid  enthalpy  hd  as, 


L{Td)=hv{Td)-hd{Td,Vv{Td)) 


evO'd )  + 


MWV 


ed(Jd )  + 


Vv(Jd) 

Pd 


(21) 


where  e  is  the  specific  internal  heat,  R0  the  universal  gas  constant,  MWV  the  molecular  weight  of 
the  vapor,  and  pv(Td)  is  the  equilibrium  vapor  pressure. 

4.1.4  Drop  Deformation 

Drop  distortion  is  important  in  the  detennination  of  breakup  in  atomization  process.  Typically 
defonnation  is  modeled  by  use  of  Taylor’s  drop  oscillator  (TAB  model)  where  distortion  is 
described  by  a  forced,  damped,  hannonic  oscillator.  The  forcing  tenn  is  given  by  the 
aerodynamic  drag,  the  damping  is  due  to  the  liquid  viscosity,  and  surface  tension  is  the  restoring 
force.  The  deformation  equation  is  written  as, 


CpPgU2  CpO  CdPi  . 

y-~r — — 5- y - — t  y 

LBPirz  PjT6  ptrz 


(22) 


where  CF,  CK,  and  Cd  are  dimensionless  numbers,  further  details  of  this  model  is  given  in  section 
4.2.4. 

4.1.5  Drop  Collisions 

The  droplet  collision  model  of  O’Rourke  and  Amsden  ( 36)  is  consistent  with  the  stochastic 
nature  of  spray  simulations,  where  only  a  subsample  of  droplets  is  tracked.  It  represents  parcels 
of  varying  number  of  drops.  The  size,  number,  and  velocity  of  parcels  detennine  the  probability 
of  colliding  with  other  parcels.  The  drop  collisions  are  accounted  via  the  source  term  fcoa. The 
probability  for  a  drop  with  index  1  to  undergo  collisions  with  a  drop  of  index  2  in  a  given 
volume  V  during  a  specified  time  interval  At  is  given  by  the  Poisson  distribution, 

Pn  —  xnexp(—x)/n\  (23) 
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where  x  =  vAt  is  the  mean.  The  collision  frequency  is  modeled  as, 

No 

v  =  y? r(ri +  r2)2||v1  -  v2||  (24) 

where  N2  represents  the  number  of  drops  in  volume  and  drop  state  2.  The  collision  in  this  model 
will  have  two  possible  outcomes,  depending  on  the  collision  impact  parameter  b;  for  cases  when 
b  <  bcr  then  the  drop  will  coalesce,  and  if  b  >  bcr  then  the  drops  undergo  an  elastic  collision 
(rebound)  maintaining  their  size  and  temperature.  The  critical  impact  parameter  bcr  —  bcr(ru  r2,  cr) 
depends  on  collision  drop  sizes  and  surface  tension,  a. 

Other  techniques  for  modeling  droplet  collision  include  the  works  proposed  by  Schmidt  and 
Rutland  (30),  it  is  based  on  the  No  Time  Counter  (NTC)  method  used  in  gas  dynamics  for  Direct 
Numerical  Simulation  calculations. 

4.2  Review  of  Liquid  Spray  Breakup  Models 

Typically  fuel  injection  is  handled  using  a  blob  injection  model  and  due  to  the  interaction  with 
the  environment  the  computational  parcels  are  subject  to  surface  instabilities. 

In  diesel  spray  applications,  the  instabilities  are  typically  described  through  KH  and  RT  models, 
which  are  used  to  predict  primary  and  secondary  breakup.  An  intact  liquid-core  breakup  length  is 
used  where  the  KH  model  alone  is  used  to  predict  primary  breakup;  downstream  of  this  critical 
length  (and  in  the  hybrid  case)  the  RT  and  KH  models  are  implemented  in  competing  manners, 
such  that  the  droplet  breaks  up  by  the  model  that  predicts  a  shorter  breakup  time.  In  the  injector 
nozzle  region,  where  droplet  velocities  are  larger,  the  RT  breakup  model  dominates,  while  in  KH 
model  is  used  further  downstream. 

The  section  will  describe  in  detail  some  of  the  leading  engine  spray  models. 

4.2.1  KH  Breakup  Model 

The  fonnulation  of  the  KH  Wave  Breakup  Model  developed  primarily  by  Reitz  and  Diwaker 
(31)  considers  a  cylindrical  liquid  jet  of  radius  a  penetrating  through  a  circular  orifice  into  a 
quiescent  incompressible  gas  chamber.  The  interaction  between  the  surrounding  gas  and  the 
liquid  jet  creates  a  number  of  infinitesimal  surface  perturbations  that  are  characterized  with 
initial  amplitude  of  17  0  and  a  spectrum  of  wavelengths  A  (k  =  2n/A),  see  figure  6.  The 
perturbation  amplitude  increases  exponentially  due  to  the  interaction  between  the  liquid  and  the 
gas;  the  growth  rate  is  complex  and  is  described  as  u>  =  <ur  +  ten,.  The  surface  perturbations  are 
described  with  the  following  relation: 

17(f)  =  R(r]0exp  [ikx  +  u>]))  (25) 
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Figure  6.  Schematic  growth  of  surface  perturbations  in  KH 
model.  Notation  1  depicts  the  liquid  phase,  while  2 
depicts  the  gas  phase  (32). 


A  dispersion  equation  relating  the  growth  rate  oo  to  the  wave-number  k  can  be  derived  assuming 
that  the  gas-phase  behaves  as  an  inviscid  fluid,  and  that  the  perturbations  are  significantly 
smaller  than  the  jet  radius  (77  «  a): 


a)2  +  2vlk2(jci 


I[(ka) 

2 kl  I-iQca)  I[(la) 

1'oika) 

k2  +  l2 10(ka )  /q (Za) 

ok 

Pi°.2 


(1  —  k2a2) 


fl2  —  k2\  I-^(kcL)  pg  1  ia)\2  2  fl2  —  k2\  I^ka^Ko^ka) 
\l2  +  k2JI0(ka)  +  pi\  k  )  \l2  +  k2)  ^(ka^K^kd) 


In  equation  26,  In  and  Kn  are  the  nth  order  modified  Bessel  functions  of  the  first  and  second 
kind,  respectively.  Primes  indicate  differentiation,  vL  is  the  kinematic  viscosity  of  the  liquid, 

U  the  gas  velocity  at  the  liquid  surface,  and  l2  —  k2  +  (jo/vt.  In  this  model  it  is  assumed  that  only 
the  fastest  growing  perturbation,  denoted  by  growth  rate  fl  and  corresponding  to  the 
wavelength  A  ,  will  have  the  strongest  impact  on  liquid  breakup.  Figure  7  shows  the  solution  of 
the  dispersion  equation  24  for  various  liquid  Weber  and  Ohnesorge  numbers.  It  relates  a  critical 
maximum  wave  growth  rate,  fl ,  to  nondimensionalized  wave  length,  A  ,  at  the  fluid  conditions 
presented  (i.e.,  a  peak  growth  rate  of  fl  =  60  at  conditions  We  —  30  ,  Oh  —  30  is  predicted  with 
a  corresponding  wavelength  of  A  =  0.5  ).  In  this  model,  fl  —  60  and  A  =  0.5  represents  the 
fastest  growing  waves  on  the  liquid  surface  responsible  for  the  breakup. 
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Figure  7.  Wave  growth  rate  vs.  wave  number  as  a  function  of 
Weber  number  and  Ohnesorge  number. 


A  simplification  is  made  by  Reitz  (32)  by  generating  curve  fits  of  numerical  solutions  to 
equation  26  from  which  the  following  expressions  were  obtained  for  the  maximum  growth  rate 
co  —  fl  and  the  corresponding  wavelength  A  —  A, 


A  (l  +  0.45Z0-5)(l  +  0.4r0-7) 

-  =  9.02  - - — - 

a  (l  +  0.8714/e]-67)  ' 


(27) 


^  (pia3\  _  0.34  +  0.38M/e]-5 
n  =  (i+z)(i  +  i.4r0-6) 


(28) 


where, 

Z  =  We?s/Rel ,  T  =  14/e]-5,  Wet  =  ptU2a/a  ,  Weg  =  pgU2a/a,  and  Ret  =  Ua/vL 


Figure  8a  and  b  show  the  behavior  of  equations  27  and  28  as  the  Weber  number  is  varied.  The 
effect  of  increasing  the  Weber  number  leads  to  increases  in  maximum  growth  rate  behavior  and 
decreases  in  wavelength.  In  contrast,  variations  in  the  Ohnesorge  number  (depicted  by  Z)  is  seen 
to  reduce  the  wave  growth  rate  and  increase  wavelength  significantly  as  the  viscosity  increases. 
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Weber  number,  We2  Weber  number,  We2 


Figure  8.  Physical  behavior  of  (a)  wavelength  and  (b)  growth  rate  of  the  most  unstable  surface  wave. 

It  is  observed  that  for  increasing  Weber  numbers  the  growth  rate  increases  and  the  corresponding 
wavelength  decreases.  This  is  consistent  with  the  experimental  observation  that  for  increasing 
injection  velocity  breakup  in  enhanced  and  the  average  diameters  of  the  resulting  droplets  are 
reduced.  Viscous  effects,  appearing  in  the  Reynolds  and  Ohnesorge  number,  serve  to  attenuate 
and  reduce  the  wave  growth  rate. 

Droplet  size  is  estimated  assuming  a  linear  relationship  between  the  resulting  droplet  radius  rd, 
and  the  wavelength  A  of  the  most  unstable  surface  disturbance  (where  B  is  of  order  unity) 


rd  =  BA 


(29) 


The  spray ’s  liquid  core  length  (considered  intact)  may  be  approximated  by  considering  the  mass 
removed  from  the  jet  by  the  atomization  process, 


L  = 


(30) 


where  f(T)  asymptotically  approaches(30,5/6)  for  (T  >  100)  which  is  typical  for  diesel  sprays. 
The  constant  c  ranges  from  approximately  15-30  and  accounts  for  various  effects  of  the  inner- 
nozzle  flow  that  are  not  resolved  in  detail. 

The  spray’s  half-angle  (a/ 2)  in  the  atomization  region  (of  high-speed  jets)  has  been  specified 
based  on  the  assumption  that  the  droplet  velocity  component  normal  to  the  spray  direction  is 
proportional  to  the  wave  growth  rate  of  the  most  unstable  wave, 
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where  the  constant  A  accounts  for  the  nozzle  geometry  and  it  has  been  defined  in  terms  of  the 
length  to  diameter  ratio  of  the  nozzle  hole, 


A  =  3.0  + 


^noz/^ 
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noz 


(32) 


4.2.2  Blob  Injection  Model 

In  diesel  engine  applications,  Reitz  (32)  and  Reitz  and  Diwaker  (33)  have  used  a  blob  injection 
model  that  continuously  injects  into  the  gas-phase  large  drops  (blobs)  with  a  diameter 
comparable  to  the  size  of  the  nozzle  hole.  The  frequency  of  the  addition  of  new  blobs  is  related 
to  the  fuel-injection  rate,  assuming  constant  density  of  the  liquid  fuel  and  ideally  spherical  blobs. 
The  KH  model  is  applied  immediately  after  the  injection  region  to  provide  the  aerodynamic 
instabilities  that  will  begin  to  grow  on  the  droplet  surface;  this  causes  smaller  secondary  droplets 
to  be  sheared  off  of  the  parent  droplet  surface  as  depicted  in  figure  9. 


Figure  9.  Illustration  of  the  blob-injection  model  of  Reitz  et  al.  (32). 

Reitz  (32)  proposed  the  calculation  of  the  resulting  droplet  radii  through  the  use  of 
equations  27  and  28)  for  the  fastest  growing  wavelength  A  and  the  maximum  growth  rate  fl,  and 
by  using  the  following  relation, 

rd  =  B0A  (33) 

where  B0  —  0.61  is  used  for  a  droplet  “stripping”  breakup  regime  characterized  as 
25  <  We  <  50.  In  cases  where  higher  injection  velocities  are  applicable,  leading  to  catastrophic 
breakup  in  the  atomization  regime  (or  We  >  50),  the  following  relation  was  proposed  by  Hwang 
et  al.  (34), 

B0  =  0.3  +  0.6P  (34) 

where  P  is  a  random  number  within  the  interval  between  zero  and  one.  In  doing  this,  a 
distribution  of  droplet  sizes  is  obtained  which  is  more  realistic  in  breakup  of  high-speed  jets.  The 
model  restricts  the  use  of  equations  27  and  28  only  for  the  case  where  B0 A  <  a,  the  resulting 
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droplet  diameter  is  less  than  the  radius  of  the  original  drop.  When  this  condition  is  not  satisfied, 
rd  is  estimated  as, 


((3na2U /2Ci)1/3 
(  (3a2A/4)1/3 


for 


B0 A  >  a,  calculated  once 


(35) 


Equation  35  is  based  on  the  assumption  that  the  jet  disturbance  has  a  frequency  of  /  =  Cl/2n 
(i.e.,  one  drop  is  formed  each  wave  period),  or  that  droplet  size  is  determined  from  the  volume  of 
liquid  contained  under  one  surface  wave. 

The  breakup  and  generation  of  new  small  droplets  lead  to  a  reduction  in  size  of  the  original  blob. 
A  mathematical  description  of  the  temporal  change  in  radius  of  the  parent  drop  is  written  as, 

da  a  —  rb 

r~  (36) 


where  r  is  the  breakup  time  given  as, 

r  =  3.726Bi  — 
1  ACL 


(37) 


and  the  constant  Bt  is  introduced  to  account  for  internal  flow  effects  of  the  nozzle  flow  on  the 
breakup  time.  Reitz  and  Diwaker  (31,  33)  suggests  using  a  value  Bx  =  20  where  other  references 
report  more  accurate  values  with  values  ranging  from  Bt  =  1.73  up  to  B±  —  30.  To  reproduce 
the  spray  cone  angle  the  child  droplets  emanating  from  the  blobs  have  a  velocity  component 
perpendicular  to  the  main  spray  orientation.  Maximum  values  of  this  component  are  obtained 
from  the  spray  half- angle  specified  in  equation  29.  It  also  suggests  using  an  even  distribution 
between  zero  and  the  maximum  normal  velocity  for  the  various  droplets  in  order  to  achieve  a 
realistic  droplet  density  within  the  spray. 


Additionally,  it  is  also  shown  that  in  the  limits  We  ->  0  and  We  ^  co  the  characteristic  breakup 
time,  r  ,  takes  the  fonn  of, 


r 


t  =  < 


v 


,  for  bag  breakup 


,  for  shear  breakup 


(38) 


Conditions  for  bag  and  shear  breakup  are  experimentally  detennined  and  reported  to  be  We  > 
6  and  We/^ReG  >  0.5  ,  respectively,  with  the  Weber  and  Reynolds  number  based  on  the 
droplet  radius. 
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4.2.3  Linearized  Instability  Sheet  Atomization  (LISA)  Breakup  Model 

Pressure  swirl  atomizers  are  often  used  in  order  to  establish  hollow  cone  sprays  as  is  typical  for 
gasoline  direct  injection  engines.  These  kinds  of  sprays  are  characterized  by  high-atomization 
efficiencies  and  effective  fuel-air  mixing  that  can  be  obtained  with  moderate  injection  pressures 
(in  the  range  of  5  to  10  MPa)  (35).  This  atomizer  accelerates  the  liquid  through  nozzles  known  as 
swirl  ports  into  a  central  swirl  chamber,  where  the  swirling  liquid  pushes  against  the  chamber 
walls  and  develops  a  hollow  air  core.  This  in  turn  produces  a  thinning  sheet  that  emerges  from 
the  orifice,  which  is  unstable,  breaking  up  into  ligaments  and  then  into  droplets.  The  process 
from  internal  injector  flow  to  a  hollow  spray  can  be  divided  into  three  steps:  film  formation, 
sheet  breakup,  and  atomization.  This  mechanism  is  clearly  depicted  in  figure  10. 


Figure  10.  Schematic  (a)  sheet  and  spray  formation  with  a  pressure  swirl  injection,  (b)  breakup  mechanism 
of  planar  liquid  sheets  (32). 


Primary  breakup  from  direct  injection  engines  utilizing  pressure  swirl  atomizers  is  typically 
modeled  using  the  LISA  technique  as  presented  by  Senecal  et  al.  (35).  The  model  emphasizes 
fluid  mechanical  principles  and  aims  at  removing  the  required  need  for  experimental  data.  The 
process  of  disintegration  of  the  liquid  sheets  into  ligaments  and  primary  droplets  is  given  in  the 
following. 


A  Bemoulli-type  expression  is  used  to  provide  a  pressure-dependent  relationship  for  the  fuel 
velocity  exiting  the  injector, 


u  =  cd 


2Apinj 

Pi 


(39) 


An  axial-velocity  component  can  be  determined  by  the  cone  half-angle  (0)  of  the  spray 
(assumed  to  be  known  for  a  given  injector): 

u  -  Ucos(0) 


A  discharge  coefficient  Cd  —  0.7  is  assumed  based  similarity  considerations  between  the  swirl 
ports  and  nozzles.  In  order  to  preserve  conservation  of  mass  from  the  measured  mass  flow  rate, 
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minj,  Reitz  (32)  suggests  to  enforce  the  following  expression  (discharge  coefficients  taking 
values  below  0.7  are  assumed  to  be  representative  of  an  air-core  existing  in  the  center  of  the 
rotating  flow) 


C, 


max 


4  mini 

0  7  \ 

Pi 

J'  ’  nd*ozpl  cos(9)  \ 

2Apinj 

(41) 


The  continuity  equation  relating  the  thickness  tf  of  the  liquid  film  inside  the  injector  to  the 
measured  mass  flow  rate  is  given  as, 

rilinj  —  TCUtf  (d-noz  tf)  (42) 

Following  the  same  assumptions  as  the  KH  wave  breakup  model,  a  spectrum  of  infinitesimal 
disturbances  imposed  on  the  liquid-sheet  surface  from  the  liquid-gas  interaction  causing  the 
amplitudes  of  the  waves  to  grow  as  in  equation  23,  p(t)  —  R(p0exp  [ikx  +  u>]).  Similarly,  the 
most  unstable  disturbance  having  the  largest  value  of  u)r,  denoted  by  H,  is  assumed  to  be 
responsible  for  sheet  breakup  and  an  expression  for  a  dispersion  relation,  <nr  =  <ur(/c),  is  clearly 
needed  to  extract  this  parameter.  Senecal  et  al.  (35)  provided  a  simplified  dispersion  relationship 
assuming  sinuous  wave  mode  disturbances  on  the  liquid  sheet  surfaces  (waves  traveling  on  both 
surfaces  of  the  liquid  sheet  travel  exactly  in  phase), 

(or  =  -2  vtk2  +  \4vfk4  +  —  U2k2  —  (43) 

4  pi  pi 

Equation  43  provides  a  dispersion  relation  (for  sinuous-wave  mode  surface  disturbances)  that  can 
be  obtained  under  the  following  assumptions:  (1)  using  an  order  of  magnitude  analysis  from  the 
inviscid  solutions  to  show  that  the  second  order  tenns  in  viscosity  can  be  neglected  in 
comparison  to  other  tenns;  (2)  for  cases  where  Weber  numbers  of  Weg  —  27/16  is  exceeded, 
short  waves  will  grow  on  the  sheet  surface,  with  a  growth  rate  independent  of  the  sheet 
thickness;  (3)  the  gas-to-liquid  density  ration  has  to  meet  the  condition,  pg/pi  «  1. 

The  breakup  time  used  to  calculate  the  critical  amplitude  of  surface  disturbances  (the  time  where 
ligaments  are  assumed  to  be  formed)  is  provided  as, 

Vb  =  V0exp(flTb) 
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where  r]b  is  the  critical  amplitude  at  breakup,  p0is  the  amplitude  of  the  initial  disturbance,  and 
H  is  the  maximum  growth  rate  that  is  obtained  by  numerically  maximizing  equation  26  as  a 
function  of  the  wave  number  k  . 

The  breakup  length  L  is  now  estimated  by  assuming  constant  velocity  for  the  liquid  sheet, 


L  =  U  rb 


(45) 


where  quantity  rib/ri0  has  been  reported  as  tjb/ri0  =  12. 

Ligament  diameter  formed  at  the  breakup  point  is  obtained  from  a  mass  balance,  assuming  that 
the  ligaments  are  formed  from  tears  in  the  sheet  once  per  wavelength.  The  resulting  diameter  is 
given  as, 


dL 


=  -A,  ■  2 h  = 


n 


(46) 


where  Ks  is  the  wavenumber  Ks  —  2n/As  corresponding  to  the  maximum  growth  rate  fi  leading 
to  the  breakup  of  the  sheet.  The  ligament  diameter  is  a  function  of  the  sheet  half-thickness  h.  at 
the  breakup  position,  which  is  related  to  its  initial  value  h0  at  the  nozzle  orifice  by, 

_  h-o  \d-noz  —  f/] 

2 L  sin(9)dnoz  —  tf  (47) 


and 


h0  *  ycos(0) 


(48) 


Further  breakup  of  ligaments  into  droplets  is  calculated  based  on  an  analogy  to  Weber’s  result 
for  growing  waves  on  cylindrical,  viscous  liquid  columns.  The  wave  number  KL  for  the  fastest 
growing  wave  on  the  ligament  is, 


3  hi 

2  /[pod 


(49) 
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Assuming  that  breakup  occurs  when  the  amplitude  of  the  most  unstable  wave  is  equal  to  the 
radius  of  the  ligament  (one  drop  will  be  formed  per  wavelength).  A  mass  balance  then  yields  an 
expression  for  the  droplet  diameter,  ddrop 


d 


drop 


(50) 


Validations  were  conducted  using  experimental  data  from  two  high-speed  automotive  fuel 
injectors  (pressure  swirl)  with  Stoddard  Solvent  as  the  working  fluid.  The  injector  is  an  inwardly 
opening  injector  with  an  air  core,  used  at  injection  pressure  —  4.86  MPa  ,  minj  =  44  mg  , 
durinj  —  3.86  ms,  and  fuel  density,  f)j  —  0.77  g /cm3 .  Figure  1  la  shows  the  qualitative  droplet 
behavior,  while  figures  lib  and  c  characterize  the  spray  penetration  and  Sauter  Mean  Diameter 
(SMD)  with  successful  comparisons  to  experiments. 


Figure  11.  (a)  Quantitative  droplet  distribution  from  atomizer;  measured  and  predicted  (b)  Spray 
penetration,  (c)  SMD  with  inwardly  opening  injectors  (Conditions:  Pinj  —  4.86  MPa  , 
minj  —  44  mg  ,  durinj  =  3.86  ms,  and  fuel  density,  pf  —  0.77  g/cm 3  ) 

4.2.4  Taylor  Analogy  Breakup  (TAB)  and  E-TAB  Breakup  Model 

O’Rourke  and  Amsden  ( 36)  have  developed  one  of  the  most  prominent  models  used  for  breakup 
calculation.  It  is  based  on  an  analogy  between  an  oscillating  and  distorting  droplet,  and  a  spring 
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mass  system.  The  restoring  force  of  the  spring  corresponds  to  the  surface  tension  forces,  external 
forces  on  the  mass  corresponds  to  the  aerodynamic  drag  forces,  and  the  damping  force  is 
represented  by  the  viscosity  of  the  liquid.  The  formulation  of  the  classical  TAB  model  is  given 
below, 

The  equation  of  the  damped,  forced  hannonic  oscillator  is  given  by: 

mx  —  F  —  ksx  —  dx  (5 1) 


where  x  is  the  displacement  (droplet  displacement  from  equilibrium),  F  are  the  external  forces 
(drag),  ks  is  the  spring  constant  (surface  tension),  and  d  is  the  damping  parameter  (liquid 
viscosity).  With  the  use  of  Taylor’s  analogy,  the  physical  dependencies  of  the  coefficients  are 
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where  CF,  CK,  and  Cd  are  dimensionless  numbers.  Additionally,  CB  is  used  to  nonnalize  x, 
y  —  x/CBr  and  normalizing  the  oscillator  equation  yields, 

Pg u2  Ckg  CdPi  . 

y  —T. - 7 - 7T - TV 

CB  Pi  r 2  ppr3  ptr2 


(53) 


with  breakup  occurring  if  and  only  if  y  >  1,  and  if  and  only  if  the  amplitude  of  oscillation  of  the 
north  and  south  poles  equals  the  drop  radius.  The  dimensionless  constants,  CF,  CK,  and  Cd,  are 
detennined  by  comparing  with  experimental  and  theoretical  results  and  have  the  following 
values:  CK  =  8,  Cd  —  0.5,  Cb  =  0.5  ,  and  CF  =  1/3.  We  also  note  that  it  has  been  reported  by 
Grover  et  al.  (37)  that  for  gasoline  sprays  a  value  of  CK  —  0.6  resulted  in  better  agreement  with 
measurements;  thus  the  need  for  model  parameter  calibration  for  different  fuels  is  clearly  needed. 


The  prediction  of  droplet  sizes  after  breakup  is  based  upon  an  energy  conservation  analysis.  In 
the  analysis,  a  balance  of  the  energy  of  the  parent  drop  with  the  energies  of  the  newly  generated 
drops  after  breakup  carried  out  to  yield, 


r  8  K  ptr3 

^2  =  1  +  20+—i> 


(54) 


where  r32  is  the  Sauter  Mean  Radius  (SMR)  of  the  parent  droplet  and  K  a  constant  that  must  be 
detennined  by  measuring  drop  sizes.  O’Rourke  and  Amsden  ( 36)  suggested  a  value  of 
K  =  10/3.  One  of  the  model  limitations  is  that  it  can  only  represent  one  oscillation  mode  (in 
reality  several  modes  exist)  corresponding  to  the  lowest  order  harmonic  whose  axis  is  aligned 
with  the  relative  velocity 
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vector  between  liquid/gas.  This  is  the  most  important  oscillation  modes,  but  for  large  Weber 
numbers  other  modes  become  more  significant  to  the  breakup  process. 

The  Enhanced-TAB  model  (E-TAB)  developed  by  Tanner  in  1997  (35)  and  Tanner  and  Weisser 
in  1998  (39)  reflects  a  cascade  of  droplet  breakup,  where  the  breakup  condition  is  calculated 
based  on  the  Taylor’s  droplet  oscillator  dynamic  analogy  (32).  The  droplet  size  is  reduced  in  a 
continuous  manner,  until  the  product  droplets  reach  a  stable  condition,  while  the  droplet 
deformation  dynamics  from  the  original  TAB  model  is  maintained. 


Breakup  occurs  when  the  droplet  distortion  (normalized  with  respect  to  the  initial  radius) 
exceeds  the  critical  value  of  1.  The  rate  of  droplet  creation  is: 


— m(t)  =  -3  Kbrm(t) 


(55) 


where  m(t)  is  the  mean  mass  of  the  droplet’s  product  distribution,  and  Kbr  a  constant  depending 
on  the  breakup  mechanism.  Further  analysis  leads  to  an  exponential  relation  between  the 
products  (r)  and  parents  (a)  droplet  radius  as, 


_  —  q  Kbut 

a 


(56) 


Tanner  (33)  TAB/E-TAB  model  performance  was  validated  through  a  range  of  conditions 
representing  low-to-intermediate  pressure  injection  into  an  ambient  nonevaporating  environment. 
The  breakup  constant  for  the  E-TAB  were  determined  through  comparison  and  experimental 
matching  of  cross-sectional  drop  size  distribution  in  the  dense  region.  Figure  12  shows 
improvements  for  droplet  velocity  and  droplet  size  in  comparison  to  the  traditional  TAB  method 
for  an  injection  time  tinj  —  1  ms  at  100  nozzle  diameters  downstream  from  injection 
(100do  =  15  mm). 
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Figure  12.  Drop  distribution  in  the  near  field  fort  inj  = 
1.0  ms  after  100  d  downstream 


( PLnj  —  300  bar). 
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Further  model  validation  is  demonstrated  with  results  in  good  agreement  with  measurements  as 
shown  in  figure  13  liquid  penetration  and  the  spray  angle  (radial  expansion).  Note  that  a  key 
feature  of  the  TAB  models  is  the  ability  to  compute  the  radial  expansion  as  part  of  the 
calculation. 


Figure  13.  Comparison  of  liquid  penetration  profiles  with  E-TAB  for  3  cases  at  Pinj  —  300  bar.  Top  case 

conditions:  Pgas  =  11  bar  ,  tinj  =  2.5  ms;  Middle  case:  Pgas  =  30  bar  ,  tinj  =  4.0  ms  ;  Bottom  case: 
Pgas  =  50  bar  ,  tinj  =  5.0  ms  . 

4.2.5  RT/Hybrid  KH-RT/Modified  KH-RT  Breakup  Model 

This  model  is  formulated  based  on  theoretical  analysis  of  the  stability  of  liquid-gas  interfaces 
when  accelerated  in  a  normal  direction  to  the  plane  (34).  Typically  RT  instabilities  tend  to 
develop  when  the  fluid  acceleration  has  an  opposite  direction  to  the  density  gradient  (the 
interface  is  stable  when  acceleration  and  density  gradients  are  aligned).  For  a  liquid  droplet 
decelerating  by  drag  forces  in  the  gas  phase  it  means  that  instabilities  may  grow  unstable  at  the 
trailing  edge  of  the  droplet, 

The  acceleration  of  a  droplet  is  due  to  the  drag  forces  and  is  mathematically  expressed  as, 


I  |  3  p n  Pyel 


(57) 


where  vrei  is  the  relative  velocity  between  the  droplet  and  gas,  r  is  the  droplet  radius.  Assuming 
a  linearized  disturbance  growth  rate  and  negligible  viscosity,  the  frequency  and  wavelength  of 
the  fastest  growing  waves  are  derived  as, 
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and, 


Art  =  2niw^7) 

It  is  observed  that  acceleration  is  the  leading  term  causing  a  rapid  growth  of  RT  instabilities, 
whereas  surface  tension  counteracts  the  breakup  mechanism. 

Traditional  (standard)  RT  models  define  the  growth  rate  as  follows, 
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More  recent  efforts  have  reformulated  the  growth  rate  to  account  for  viscosity  variations  (14), 
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The  breakup  time  is  found  as  the  reciprocal  of  the  frequency  of  the  fastest  growing  wave, 

‘-=n  (62) 

The  size  of  the  product  droplet  is  calculated  based  on  the  RT  wavelength  A  ,  and  the  breakup 
process  is  only  allowed  when  A  <  dparent.  Typically  the  RT  model  is  used  in  combination  of  the 
KH  model  as  a  hybrid  KH-RT  model  to  describe  secondary  droplet  breakup.  In  the  hybrid  case, 
the  RT  and  KH  models  are  implemented  in  competing  manners,  such  that  the  droplet  breaks  up 
by  the  model  that  predicts  a  shorter  breakup  time  (40).  In  the  injector  nozzle  region,  where 
droplet  velocities  are  larger,  the  RT  breakup  model  dominates,  while  in  KH  model  is  used  further 
downstream. 

An  additional  constraint  is  enforced,  that  is  used  in  order  to  reproduce  experimentally  intact  core 
or  breakup  lengths, 


the  RT  breakup  model  would  predict  very  rapid  breakup  directly  at  the  nozzle  exit;  thus  it  is 
switched  off  within  the  breakup  length  provided.  Since  the  RT  model  predicts  the  disintegration 
of  a  parent  droplet  into  a  number  of  equally  sized  product  droplets,  the  combination  of  RT-KH 
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model  counteracts  the  formation  of  sprays  with  bimodal  droplet  size  distributions  as  they  would 
be  predicted  if  the  KH  model  alone  is  applied  as  the  sole  mechanism. 

The  modified  KH-RT  model  removes  the  restriction  of  enforcing  a  liquid  breakup  length  “a 
priori.”  Primary  breakup  (of  the  injected  drops)  is  treated  with  KH  instability  model,  while 
secondary  breakup  is  modeled  by  competing  effects  of  the  KH  and  RT  mechanisms  (14). 

Figure  14  shows  validations  conducted  with  the  Hybrid  KH-RT  model  for  the  nonevaporating 
case  at  two  rail  pressures  of  Praii  —  170,  210  bar.  Good  agreement  is  obtained  with  liquid 
penetration  profiles  particularly  at  the  liquid  dense  region  while  using  a  Caterpillar  HEUI  3 15B 
injection  system  with  Viscor  and  cerium  blend  as  injection  liquid  fuel.  On  figure  14a,  note  that 
higher  rail  pressures  shows  lower  penetration  rates,  and  hence  lower  spray  penetration,  for  both 
experiments  and  simulations.  Also  note,  on  figure  14b,  that  more  accurate  rate  of  injection  (ROI) 
measurements  are  obtained  with  x-ray  measurement.  Comparisons  with  x-ray  ROI  and  Bosch 
rate  meter  ROI  show  gross  underpredictions  in  liquid  penetration  with  the  ROI  meter. 


Figure  14.  KH-RT  liquid  model  validation  (a)  penetration  length  validation  with  x-ray  data  (b)  effect  of  ROI  on 
penetration  rates. 

It  is  noted  that  different  time  breakup  constants  are  used  depending  on  test  conditions 
B1  —  20  (evaporating)  and  —  40  (nonevaporating),  while  the  size  constant  remains  as 
B0  —  0.61.  Similarly  the  RT  constants  will  vary  as  CRT  —  0.1  (evaporating)  and 
CRT  —  0.2  (nonevaporating);  Cx  —  10  (evaporating)  and  Cx  —  20  (nonevaporating);  with  RT 
time  constant  fixed  at  CT  —  1.0  . 

4.2.6  KH-ACT  (Aerodynamics,  Cavitation,  and  Turbulence)  Breakup  Model 

In  2010,  Som  and  Aggarwal  (41)  introduced  a  primary  breakup  model  that  accounts  for  the 
effects  of  cavitation,  turbulence  and  aerodynamics  arising  from  the  injector  nozzle  flow.  In  this 
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study,  the  effects  of  the  internal  turbulent  channel  flow  emanating  from  the  nozzle  were  reported 
to  be  important  in  predicting  a  more  accurate  rate  of  atomization,  liquid  penetration  depth,  and 
affecting  the  size  of  the  child  droplets. 

The  model,  KH-ACT  is  quasi-dynamically  linked  to  the  inner  nozzle  flow  simulation  to  provide 
the  cavitation  and  turbulence  transient  levels  at  the  nozzle  exit  as  initial  and  boundary  condition 
conditions  as  follows. 

The  turbulence  induced  breakup  model  consists  of  calculating  the  turbulent  time  and  length  scale 
assuming  isotropic  internal  flow  and  using  the  classical  decay  laws, 
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(64) 


where  K  (t)  and  e(t)  are  the  instantaneous  kinetic  energy  and  turbulent  dissipation  rate,  and 
6), and  Ce  are  the  model  constants.  Based  on  an  isotropic  turbulence  assumption,  and  neglecting 
the  diffusion,  convection,  and  production  terms  in  the  k  —  e  equation  the  following  estimates  are 
used  for  K(t)  and  e(t), 
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(66) 


where  K0  and  e0  are  the  initial  values  at  the  nozzle  exit  at  start  of  injection  (SOI),  detennined 
from  precursor  nozzle  flow  simulations  conducted  by  Som  and  Aggarwal  (41). 


The  cavitation  induced  breakup  model  considers  the  generation  of  cavities  in  the  nozzle  injector 
that  may  travel  downstream  and  reach  the  nozzle  exit;  the  mechanism  of  implosion  enhances  jet 
atomization  process  decreasing  characteristic  breakup  timescales.  The  timescale  is  set  by  the 
faster  of  two  physical  mechanisms  occurring  from  the  burst  of  the  cavity  at  the  jet  periphery  or 
by  collapse  before  reaching  it, 


T-CAV  m^n{j collapse'- T burst) 


(67) 


The  bubble  collapse  time  is  calculated  from  the  Rayleigh  Plesset  theory  as  the  time  taken  for  a 
bubble  of  a  given  radius  “r”  to  decrease  to  zero, 
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T- collapse 


0.9145  Rcav 


(68) 


where  pv  is  the  fuel  vapor  pressure,  pt  is  the  fuel  density,  and  RCAV  is  the  effective  radius  of  an 
equivalent  bubble  from  the  nozzle,  RCAV  —  rhoiey]  1  —  Ca;  and  the  area  reduction  coefficient 
Ca  is  calculated  from  precursor  internal  flow  nozzle  calculations,  and  rfto/e  is  the  exit  radius. 


The  average  time  required  for  a  cavitation  bubble  to  reach  the  periphery  of  the  jet  can  be 
estimated  as, 


_  rhoie  Rcav 
^  burst  —  ~l 
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(69) 


Where  the  turbulent  velocity  u'turb  —  yj2K{t)/3  is  obtained  from  the  precursor  calculation;  the 
length  scale  for  the  cavitation  induced  breakup  is  calculated  as  LCAV  —  RCAV  . 


Aerodynamically  induced  breakup  model  KH  model  is  used  to  calculate  the  instantaneous  length 
and  timescales  for  every  parcel  as, 


LkH  ~  r  rKH 

(70) 

3.2765^ 

(71) 

Q-KhA  KH 

The  ratio  of  length  and  timescales  for  each  process  is  calculated  since  the  largest  value  will 
determine  the  dominant  breakup  process, 


La  _  {LKH(t)  LCAV(t)  LT(t )~j 
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(72) 


In  the  case  where  aerodynamic  induced  breakup  is  dominant,  then  the  previously  presented  KH 
model  is  employed  for  primary  atomization.  When  cavitation  or  turbulence  are  the  dominant 
mechanisms  the  following  breakup  law  is  used, 


dr  La 
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(73) 


where  CrcAV  is  a  model  constant. 
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Figure  15  shows  the  effect  of  each  physical  process  in  the  context  of  a  nonevaporating  spray  at 
diesel  operating  conditions.  Figure  15a  represents  the  transient  behavior  of  the  liquid-penetration 
profiles  and  clearly  showing  that  injector- flow  turbulence  plays  a  major  role  in  the  atomization 
process.  Figure  15b  shows  similar  behavior  with  the  SMD  profiles,  although  the  variations  are 
not  as  pronounced  as  in  the  penetration  profiles. 


Time  (ms) 


Figure  15.  Primary  breakup  model  effect  on  nonevaporating  (a)  liquid  penetration  profdes  and  (b)  SMD. 

Conditions:  Pinj  =  1100  bar  ,  Pf  —  865.4  kg/m3  ,dinj  =  169  pm,  px  —  34.13,  Tm  —  298  K, 
10%0  JV2. 


Validation  with  evaporating  sprays  was  also  conducted  and  is  presented  in  figure  16.  The 
experimental  conditions  correspond  to  a  constant  volume  combustor  and  data  was  reported  for  a 
range  of  parameters  including  ambient  gas  densities  and  chamber  operating  conditions.  It  is 
observed  that  the  overall  behavior  is  captured  well  with  both  models.  However,  we  can  see 
improvements,  due  to  the  enhanced  atomization  rate  accounted  for  in  the  KH-ACT  model  that 
decreases  the  liquid  penetration  and  droplet  size. 
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Figure  16.  Comparison  between  hybrid  KH  model  and  KH-ACT  model  for  (a)  liquid  penetration  length  variations 
with  chamber  temperature  and  (b)  vapor  penetration  length  variations  with  time.  Conditions  for  the 
orifice  diameter,  injection  pressure,  and  fuel  temperature  are:  246  /im,  1420  bar,  and  438  K, 
respectively. 


5.  Summary  and  Outlook 


This  report  has  covered  important  modeling  aspects  in  diesel-engine  sprays.  In  this  approach, 
one  models  the  breakup  behavior  and  solves  the  fundamental  transport  equations  that  govern  the 
evolution  of  the  one-point,  one-time  probability  density  function  for  a  set  of  variables  that 
determines  the  state  of  the  process.  The  emphasis  has  been  on  atomization  models  that  are  often 
used  with  Lagrangian  particle/Eularian  numerical  solution  solvers  (3D-CFD).  Important 
developments  relative  to  this  baseline  have  been  discussed;  this  also  includes  a  review  of  the 
relevant  literature,  and  the  physical  submodels  for  calculating  droplet-coupled  interaction  with  its 
surroundings  (i.e.,  acceleration,  radius,  temperature,  deformation,  and  collision),  this  is 
summarized  in  table  1 . 
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Table  1.  Physical  droplet  models  used  with  Lagrangian  calculations. 


Droplet  Model 

Mathematical  Description 

Equation  Number 

Acceleration 

<1  ^  Og  ||Vr|| 

dty  =  8%,  r  V'  +  g 

Equation  13 

Radius 

^-r2  =^DvBdShd 

cit  Pd 

Equation  14 

Temperature 

dT ri  diTifj 

Cdmd  ~  Ph^d  +  L(Td) 

Equation  18 

Deformation 

Cf  Pg  n2  C Ka  CdPi  . 

y  =  — - 7 - ry - —  y 

CB  Pi  rz  ppr6  ptrA 

Equation  22 

5.1  Summary 

Significant  success  has  been  achieved  in  the  spray  modeling  realm  in  terms  of  capturing  the 
breakup  process  in  reacting  and  nonreacting  systems.  The  number  of  publications  and  validated 
spray  models  is  a  strong  indicator  of  its  success  having  a  wide  range  of  applicability.  Indeed, 
these  models  have  played  a  strong  role  in  enabling  new  technologies  and  optimizing  engine 
configurations.  However,  the  literature  also  reveals  a  significant  amount  of  empiricism  in 
atomization  models  that  often  leads  to  model  calibration  (i.e.,  tuning)  when  conducting 
validation.  Variations  in  physical  properties  of  the  model  (change  of  fuel),  boundary  conditions 
(injection  properties),  and  changes  from  low-to-high  fidelity  models  have  been  shown  to  require 
additional  recalibration  (18, 19).  Models  tend  to  be  developed  with  numerical  constants  that 
allow  control  of  a  specific  physical  process,  such  as  penetration  length,  droplet  size  (SMD), 
breakup  length,  etc.,  hence  promoting  calibration.  Surprisingly,  diesel  engine  spray  models,  such 
as  the  ones  reviewed  in  this  article,  often  lack  a  complete  description  of  its  physical  mechanisms 
critical  to  its  behavior.  While  cavitation  and  turbulence  effects  play  an  important  role  in  the 
atomization  of  high-pressure  fuel  injectors,  breakup  models  have  not  reached  a  mature  state  in 
tenns  of  its  description.  Its  implementation  in  a  contemporary  spray-atomization  model  has 
always  presented  several  technical  challenges  and  has  not  been  fully  achieved  to  date.  Table  2 
illustrates  some  of  the  important  features  of  the  atomization  models  presented. 
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Table  2.  Liquid  spray  atomization  modeling  table. 


Atomization 

Models 

Representative  Model  Equations 

Model  Constants 

Model  Features 

KH  (31) 

tkh  ~  3.726B1  —  ,  rd  =  B0A KH 

B0  =  0.61 

1.73  <  B1  <  30 

Primary’  Wave  Model 
(aerodynamic  effects) 

Blob  Injection  (33) 

~(-)a 

Rosin-Rammler:  f(r)  —  af~ara~1e  w 

Chi-Squared:  f(r)  —  f1e  (?) 

Rosin-Rammler,  x 2 
droplet  distributions 

Primary’  Wave  Model  (with 
statistical  droplet  distribution) 

LISA  (35) 

a>r  =  — 2v;fc2 

rb  =  iln 

+ 

A 

IVbj 

Lvffc4  +—U2k2  -  — 

Pi  Pi 

\  rl 

)  i  “drop  —  1^“ J 

No  constants, 
based  on  fluid 
mechanical  principles 

Sheet  Atomization  Model 
(pressure  swirl  applications) 

TAB  (36) 

..  CFPgU2  CKa  Cdm  . 
y  =  2  j  2  y 

CB  p,  r2  p,r3  p,r2 

r  8  K  p;r3  /6  K  —  5\ 

r32“1+  20  +  a  y(  120  J 

CK  —  8,  Cd  =  0.5, 

Cb  —  0.5  ,  and 

CF  =  1/3 

Primary’  Droplet  Distortion 
(with  WA  VE  model) 

E-TAB  (38,39) 

r/a  =  e~Kbut 

Kbu 

Primary’  Droplet  Distortion 
Model  (with  WA  VE  effects) 

RT  (14,3 1,40) 

flRT  — 

Ur\\\F\(Pi-Pa)]1/4 

3  3ct 

Depends  on 

WAVE  model 

Primary/Secondary  fluid 
acceleration  model 

KH-RT  (40) 

0)RT  — 

'  (Pl-Pg\  klr° 

k RT  1  )  CL 

\Pl  +  P#/  Pi  +  P# 

Depends  on 

WAVE  model 

Primary’  Hybrid  Model 
( WA  VE/Acceleration) 

Modified  KH-RT 

2  (Pl+Pg\ 

kRT  | 

'Pi  ~  Pg  \  „  ,  ,.4  (Vi+Ug'f 

Depends  on 

WAVE  model 

Primary  Hybrid  Model 
( WA  VE/Acceleration) 
**Breakup  length  is  not  enforced 

U)RT  —  KRT  1  T 

\Pl+PgJ  y 

1  Ct  1  E  nT  1  1 

,Pz+Pfl/  Pz+Ps  \Pz  +  PgJ 

KH-ACT  (41) 

„  _  rhole-RCAV  _  _  3.276 Bjr 

T burst  T KH  n  . 

^0  ^  Uturb  CLKHA KH 

K(t)>  e(0<  PcAV 

1.73  <  B1  <  30 

Primary’  Breakup  Model 
( WA  VE/Turbulence/Cavitation) 

5.2  Unresolved  Issues  and  Directions 

Cavitation  develops  due  to  pockets  of  low-static  pressure  and  abrupt  changes  in  the  internal 
nozzle  geometry.  Turbulence  is  generated  from  vapor  phase  disturbances,  internal  flow 
conditions,  wall  roughness,  and  any  microscale  manufacturing  imperfections  on  the  surface. 
Cavitation  will  increase  the  liquid  velocity  at  the  nozzle  exit,  and  hence  will  influence  the 
fonnation  and  quality  of  the  emerging  spray.  Improved  spray  development  is  believed  to  lead  to 
a  more  complete  combustion  process,  lower  fuel  consumption,  and  reduced  emissions.  However, 
cavitation  can  decrease  the  flow  efficiency  (discharge  coefficient)  due  to  its  effect  on  the  exiting 
jet.  Also,  imploding  cavitation  bubbles  inside  the  orifice  can  cause  material  erosion,  thus 
decreasing  the  life  and  perfonnance  of  the  injector  (29).  This  effect  has  only  been  directly 
accounted  for  in  the  KH-ACT  model  developed  by  Som  and  Aggarwal  (41)  and  introduced  in 
previous  sections.  Here  an  empirical  approach  was  utilized  to  account  for  turbulence  and 
cavitation,  based  on  idealized  assumptions  such  as  internal  isotropic  turbulent  flow  and 
approximations  in  cavitation  number  receiving  some  success  for  reacting  sprays. 

Several  numerical  studies  exist  with  the  intent  of  modeling  the  flow  inside  injector  under  diesel 
operating  conditions  (42-45).  Typically  several  multiphase  techniques  such  as  interphase 
tracking  methods,  two  fluid  nozzle  flow  models,  and  continuum  flow  models  or  homogenous 
equilibrium  models  (HEMs)  have  been  utilized  to  study  the  effects  of  cavitation  in  diesel  injector 
nozzle.  Recently  in  2013,  Payri  et  al.  (42)  has  achieved  significant  success  in  predicting  injector 
quantities  for  a  six  orifice  microsac  nozzle  with  cylindrical  holes,  and  inclined  to  cavitation 
effects  using  continuum  flow  models.  A  similar  model  was  used  by  Salvador  et  al.  (45)  with 
high-fidelity  LES  techniques  to  model  cavitation  in  multihole  microsac  nozzles.  Validation  was 
conducted  in  terms  of  mass  flow,  momentum  flux,  and  effective  velocity  under  real  operating 
conditions  with  success  and  the  mutual  interaction  between  cavitation  and  turbulence  presented. 
Flow  turbulence  and  vorticity  is  reported  to  affect  the  appearance  of  cavitation  near  the  orifice 
inlet,  and  once  vapor  phase  is  present  it  will  promote  the  turbulence  production  since  the  highest 
values  of  vorticity  in  the  nozzle  are  found  at  the  interphase  (45).  Several  other  numerical  studies 
exist  using  Lagrangian  multidimensional  CFD  models  to  account  for  the  onset  and  development 
of  cavitation  where  good  agreement  has  been  reported  including  higher  order  diagnostics  such  as 
turbulent  velocity,  which  is  equally  influenced  by  the  inherent  unsteadiness  from  cavitation.  It 
has  been  concluded  that  cavitation  modeling  has  reached  a  stage  of  maturity  at  which  it  can 
consistently  identify  many  of  the  effects  of  nozzle  design  on  cavitation,  thus  making  a  significant 
contribution  to  nozzle  performance  and  optimization  (46,  47). 

Recently,  several  studies  invoking  the  use  of  DNS  methodology  have  surfaced  as  a  tool  to 
develop  next-generation  primary  models  (48-53).  Particular  attention  is  given  to  the  lack  of 
theories  relating  to  the  detennination  of  the  smallest  length  scale  in  turbulent  liquid  gas  flows.  In 
several  works,  the  DNS  methodology  has  been  utilized  to  provide  new  databases  and  promote 
the  development  of  an  LES  formulation  of  liquid  atomization  (53)  through  the  use  of  the  ELSA 
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model.  In  ELSA,  a  single-phase  Eularian  approach  is  used  to  describe  the  primary  breakup 
region  (54).  In  the  downstream  region,  where  the  spray  is  considered  to  be  diluted  enough,  the 
spray  is  represented  through  conventional  Lagrangian  particles  and  their  trajectories  and 
evolutions  are  tracked.  It  was  reported  that  this  approach  needs  only  two  modeling  constant,  and 
comprehensive  validation  is  found  for  the  whole  set  of  operating  conditions  without 
manipulating  the  constant  values  (51,  54).  Similar  efforts  were  carried  out  by  Desjardin  et  al. 
(50),  where  he  modeled  a  small  portion  of  the  breakup  region  of  a  simplified  diesel  injector. 
Through  the  use  of  advanced  numerical  techniques,  Desjardin  was  able  to  show  the  robustness 
and  accuracy  of  his  approach  in  modeling  atomization.  Ning  and  Reitz  (51)  extended  the 
operating  conditions  and  conducted  a  validation  study  using  ELSA  of  nonvaporizing  diesel 
sprays  achieving  successful  comparison  with  standard  KH-RT  breakup  models  and  experimental 
measurements. 

In  light  of  these  advancements,  atomization  models  have  the  opportunity  to  enhance  its  modeling 
features  to  address  its  shortcomings.  Areas  where  spray  breakup  models  are  expected  and 
recommended  to  improve  are  as  follows: 

•  Its  ability  to  accurately  describe  the  entire  spray  process  including  the  primary  breakup 
“dense”  region  to  the  “dilute”  region. 

•  Its  ability  to  capture  the  effects  of  injector  nozzle  turbulence  and  its  mutual  interaction  with 
cavitation  phenomena. 

•  Its  ability  to  resolve  the  flow-field  with  high  resolution  as  in  LES/DNS  models,  without  a- 
priori  requirements  of  recalibration. 

•  Its  ability  to  capture  variations  in  fuel  physical  properties,  without  a-priori  requirements  of 
recalibration,  important  in  studying  fuel  effects  on  atomization. 

•  Its  ability  to  include  the  effects  of  needle  deformation,  as  it  has  been  reported  that  for  high- 
injection  pressures,  needle  deformations  can  be  of  the  same  order  than  needle  lifts. 

•  Its  ability  to  describe  the  effects  of  an  electrical  field  or  charge,  important  for  engineering 
level  electrostatic  spray  applications. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AMR 

adaptive-mesh-refinement 

ARL 

U.S.  Army  Research  Laboratory 

CERFACS 

Centre  Europeen  de  Recherche  et  de  Formation  Avancee  en  Calcul  Scientifique 

CFD 

computational  fluid  dynamics 

CPF 

Constant-Pressure  Flow 

CVP 

Constant- Volume  Preburn 

DNS 

Direct  Numerical  Simulation 

ECN 

Engine  Combustion  Network 

EGR 

exhaust-gas  recirculation 

ELSA 

Eularian  Lagrangian  Spray  Atomization 

E-TAB 

Enhanced-TAB  model 

HEM 

homogenous  equilibrium  models 

IFPEN 

French  Petroleum  Institute 

KH 

Kelvin  Helmholtz 

LES 

Large-eddy  simulation 

LISA 

Linearized  Instability  Sheet  Atomization 

MTU 

Michigan  Technological  University 

NTC 

No  Time  Counter 

RANS 

Reynolds  Average  Navier  Stokes 

ROI 

rate  of  injection 

RT 

Rayleigh  Taylor 

SMD 

Sauter  Mean  Diameter 

SMR 

Sauter  Mean  Radius 

SNL 

Sandia  National  Laboratories 
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SOI 

TAB 

TU/e 


start  of  injection 
Taylor  Analogy  Breakup 
Eindhoven  University  of  Technology 
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